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This  report  describes  a  computer  program  which  may  be  used  to  in- 
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vestigate  the  characteristics  of  a  multicavity  klystron  ’  ’  ,  It  is  written  for 

the  IBM  7090  in  "Fortran"  language  so  that  it  may  be  intelligible  to  people 
other  than  the  originator  of  the  program.  It  is  written  for  solid  or  hollow 
beam  klystrons.  For  maximum  flexibility  it  is  in  the  form  of  several  subrou¬ 
tines  connected  by  a  master  control  program.  It  is  complete  in  the  sense  that 
it  leaves  nothing  for  the  user  to  compute  manually. 
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I.  Introduction 


The  following  describes  a  multicavity  klystron  (MCK)  computer  pro¬ 
gram  which  has  been  written  for  the  IBM  7090  computer  at  Lincoln  Laboratory. 
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It  is  based  upon  two  papers  ’  describing  an  improved  MCK  analysis  and  a 
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third  in  which  a  potential  depression  correction  is  derived. 

It  has  been  assumed  that  the  reader  is  familiar  with  the  above  mentioned 
reports,  and  this  manual  is  intended  only  as  an  aid  to  the  use  of  the  program. 

The  program  was  written  in  the  IBM  709  "FORTRAN"  language  which, 
in  addition  to  being  a  compact  and  efficient  language,  provides  an  excellent 
record  of  the  computation  which  is  easily  understood  by  people  other  than  the 
originator  of  the  program.  The  writers  have  sought  to  achieve  three  qualities 
in  the  program:  generality,  adaptability  to  future  changes,  and  completeness. 

The  program  achieves  generality  insofar  as  it  is  good  for  any  solid  or  hollow 
beam  klystron  with  from  1  to  6  {this  number  can  be  extended  by  altering  input 
and  output  format  statements)  cavities  of  equal  or  unequal  spacing. 

The  program  is  flexible  and  adaptable  to  future  changes  since  it  has 
been  written  in  the  form  of  several  subroutines  connected  by  a  master  control 
program.  Any  future  change  will  merely  require  altering  or  adding  subrou- 
tines  without  disturbing  the  rest  of  the  program. 

The  program  is  complete  in  the  sense  that  it  leaves  nothing  for  the  user 
to  compute  manually.  This  can  be  seen  from  the  following  outline  in  which  only 
the  most  elementary  parameters  of  the  tube  are  required  as  input  to  the  program. 


’i'No  attempt  has  been  made  to  optimize  the  program  for  the  sake  of  brevity  or 
economy  of  operation,  nor  is  such  optimization  planned.  The  writers  realize 
that  several  correctable  redundancies  in  form  do  exist,  and  that  some  ineffi¬ 
cient  methods  are  employed.  Users  of  the  program  may  find  it  to  their  ad¬ 
vantage  to  make  several  minor  alterations. 


II.  Input 


The  input  data  are  of  two  kinds.  First,  there  are  those  which  put  into 
storage  the  mechanical  and  electrical  tube  parameters  and  associated  quanti¬ 
ties  necessary  to  the  computation.  Table  I  lists  and  defines  these  parameters. 
Secondly,  there  are  those  inputs,  listed  in  Table  II,  which  determine  the  qual¬ 
ity  and  quantity  of  the  output  data.  Some  discussion  of  these  is  necessary. 

The  quality  of  the  output  is  determined  to  some  degree  by  the  number 
of  modes  considered,  although  it  has  been  found  for  typical  tubes  that  the  effect 
due  to  modes  of  order  higher  than  the  fourth  contribute  little  to  the  end  results. 
The  decision  of  the  degree  of  accuracy  desired  has  been  left  to  the  user  in  that 
the  number  of  modes  to  be  used  in  the  calculation  is  an  input  parameter. 

Since  the  output  is  generally  in  the  form  of  some  parameter  as  a  func¬ 
tion  of  frequency,  the  number  of  frequencies  at  which  data  is  desired  partially 
determines  the  amount  of  output  data.  Output  data  may  be  obtained  over  a  16% 
or  smaller  bandwidth  for  as  many  frequencies  as  are  desired.  The  input  nec¬ 
essary  to  determine  the  bandwidth  and  number  of  frequency  points  is;  the  frac¬ 
tional  deviation  from  the  center  frequency  of  the  initial  (lowest)  operating  fre¬ 
quency  at  which  data  is  desired,  the  fractional  frequency  increment,  and  the 
total  number  of  points  desired. 

There  are,  in  addition,  three  other  parameters  not  directly  used  in  the 
computation.  These  will  have  value  1  or  0  and  will  determine  which  of  several 
options  in  output  are  to  be  chosen. 

It  might  be  advantageous  for  a  particular  tube  to  change  the  values  of 

•y^  at  which  the  searches  for  solutions  start  and  the  spacing  of  the  points  used 

~  221/2 

in  the  search.  (Note  \^  =  —  k  )  '  ,  where  are  the  roots  of  the  propaga¬ 

tion  equation  and  k  =  w/c.  )  If  the  range  is  too  wide  or  the  initial  spacing  too 
fine,  computer  time  is  wasted,  if  the  range  is  too  narrow  or  the  initial  spacing 
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too  coarse,  a  root  may  be  missed.  As  the  program  is  presently  set  up,  the 
search  for  starts  at  0.  6/3^  and  increases,  and  the  search  for  y  starts  at 
1.  8/3^  and  decreases.  This  change  may  be  effected  by  altering  statement  num¬ 
ber  121  and  the  immediately  following  (unnumbered)  statement  which  reads 
presently: 

121  GPOS(I)  =  0.6«BE(1) 

GM1NUS{I)  =  1.  8*BE(1) 

(Note  /3g  =  (o/v.  )  This  is  adequate  for  conventional  multicavity  klystrons  where 

v/c  is  of  the  order  of  0.  5  and  where  the  are  spaced  on  the  order  of  +  0. 1 

from  B  (i.e.,  w  /w  ='0.1).  Similarly,  the  initial  spacing  of  the  points  used  to 
e  q 

search  for  the  root  can  be  changed  by  altering  statement  number  153  which  reads 
presently: 

153  DDG  =  .  001 

Increasing  the  value  of  "DDG"  widens  the  initial  spacing. 
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INPUT 


Fortran  Symbol 
NCVS 
VO 
FIO 
FO 
A 
B 
C 

DGAP 

RQ(I) 

QO(I) 

QEXT(I) 

S(I) 

THETA 

FNU(I) 


TABLE  I 

TUBE  PARAMETERS 


Parameter 


a 

b 


(R/Q) 

(i) 


(i) 


Q 

Q 


(i) 


ext 


f  (i) 


ext 

.(i) 


Definition 

Number  of  cavities 

DC  voltage 

DC  current 

Center  frequency 

Drift  tube  inside  radius  (inches) 

’i'Beam  outside  radius  (inches) 

*Beam  inside  radius  (inches) 

Cap  half  length  (inches) 

(R/Q)  value  for  each  cavity 

Circuit  loading  for  each  cavity 

External  loading  for  each 
cavity 

Center  gap  to  center  gap  dis¬ 
tance  for  each  gap  pair 

Phase  shift  in  external  circuits 

Fractional  detuning  for  each 
cavity 


*In  this  treatment,  infinite  magnetic  field  is  assumed,  so  that  there  is 
no  variation  of  beam  dimensions. 


4 


TABLE  II 


OUTPUT  OPTIONS 


Fortran  Symbol 


Parameter 


Definition 


N 

SFRIQ 


FFDEL 


f 

o 


Number  of  modes 


Fractional  deviation  of  input 
frequency  (f)  from  center 
frequency  (f^)  (SFRIQ  >  -  ,  08) 

Fractional  frequency  increment 


NKK 


POTD 

SKIP 


WOT 


Number  of  frequencies  at  which 
information  is  desired 
(Note:  SFRIQ  +  NKK  «  FFDEL 
must  be  less  than  or  equal 
to  +  .  08) 

=  0  do  potential  depression 
/  0  omit  potential  depression 

=  0  do  complete  computation 
thru  gain  and  phase  as  a 
function  of  frequency 
/  0  do  all  calculations  except 
gain  and  phase  as  a  function 
of  frequency 

4  0  print  all  output  parameters 
=  0  print  only  gain  and  phase 
information 
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III.  Output 


The  heading  or  the  first  nineteen  lines  of  the  printed  output  which  are 
always  printed  verify  some  of  the  input  parameters  and  serve  to  identify  the 
tube  and  conditions  being  examined.  Each  of  these  lines  has  a  label.  The 
contents  of  these  lines  are: 

1.  The  contents  of  identification  input  card  number  1 

2.  Number  of  cavities 

3.  Number  of  modes 

4.  DC  voltage 

5.  DC  current 

6.  Center  frequency 

7.  A  - 

8.  B  - 

9.  C  - 

Lines  10  through  16  give  the  following  parameters  for  each  of  the 
cavities  considered  and  the  drift  tube. 

10.  Cavity  numbers 

11.  DC  cavity  voltage  (depressed  potential) 

12.  D 

13.  R/Q 

14.  QO 

15.  QEXT 

16.  S 


The  information  for  the  beam  in  the  drift  tube  is  presented  with  the  cavity 
data  in  the  last  column  (i.  e.  ,  as  cavity  number  NCVS  +  1). 
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The  heading  is  completed  by  printing  out  the  three  input  parameters: 

17.  NKK 

18.  FREQ  =  (SFRIQ) 

19.  DEL  =  (FFDEL) 

The  second  portion  of  the  output  is  controlled  by  the  WOT  input  option 
parameter.  It  will  be  printed  out  if  WOT  4  consists  of  major  groupings 

for  each  frequency  and  minor  groups  for  each  mode  at  each  particular  fre¬ 
quency.  The  major  group  is  denoted  by  the  printed  line: 

"FDEL  =  "xyz.  ..."  where  "xyz"  is  (frequency  -  This  will 

be  followed  by  as  many  groups  of  the  following  twenty -eight  parameters,  as 
there  are  modes  to  be  considered.  Each  minor  group  is  headed  by  the  line 
"M  =  x"  where  x  is  the  mode  number.  Each  parameter  is  given  in  tabular 
form  for  each  cavity.  The  parameters  again  have  a  label.  They  are: 


Fortran  Symbol 

Parameter 

Definition 

BE 

/I 

B  + 

Fast  wave  longitudinal  plasma 
wave  number 

B- 

P-n 

Slow  wave  longitudinal  plasma 
wave  number 

G+ 

^+n 

Fast  wave  transverse  E.M. 
wave  number 

G- 

^-n 

Slow  wave  transverse  E.M. 
wave  number 

T-f- 

+n 

Fast  wave  transverse  plasma 
wave  number 

T- 

T 

-n 

Slow  wave  transverse  plasma 
wave  number 

BQPOS 

^q+n 

BQMIN 

^q-n 
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Fortran  Symbol 

Parameter 

WQ+/WP 

o  ,  /to 
q+n'  p 

WQ-/WP 

to  /to 

q-n'  p 

WQ+/W 

to  ,  /to 
q+n' 

WQ-/W 

c 

1 

cr 

3 

ZETA  + 

ZETA  - 

Y(p+/G(p 

Y0  -  /G0 

u 

o 

c 

1 

o 

>- 

E+ 

E+n 

E- 

E 

-n 

F+ 

F- 

F 

-n 

C+ 

C+n 

C- 

C_n 

M+ 

M, 

+n 

M- 

M 

-n 

H+ 

+n 

H- 

H 

-n 

GELGO 
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At  the  end  of  each  frequency  group,  there  are  printed  out  (if  WOT  ^  0) 
the  following: 


Fortran  Symbol 

Parameter 

Definition 

GELGOP 

Gel/°o=I 

n 

(G  ,/G  ) 
el'  o'n 

QEL 

Qel 

Q  of  electronic  load 

Q 

Q 

total  Q 

QOEL 

^oel 

Q  of  cavity  and  electronic  load 

GEL 

^el 

electronic  conductance 

REL 

\l 

1/G,i 

if  WOT  /  0  and  SKIP  = 

C 

Y 

Q(i)/Q0) 

ZGR 

Re  (Z^^h 
g 

Real  part  of  gap  impedance 

ZGI 

Im(Z^*^) 

g 

Imaginary  part  of  gap  imped¬ 
ance 

FIVR 

Real  part  of  over-all  trans¬ 
admittance 

FIVI 

Imaginary  part  of  over-all 
transadmittance 

GVR 

Re  (G^) 

Real  part  of  over-all  voltage 
gain 

GVI 

Im  (G^) 

Imaginary  part  of  over-all 
voltage  gain 

After  all  of  the 

major  (frequency) 

groupings  are  printed  out,  the  cavity 

detunings  are  given.  These  are  labeled  NU  and  given  for  each  cavity.  Finally, 

the  numerical  gain,  gain  in  db,  and  phase  are  given  as  a  function  of  frequency. 
These  quantities  are  unlabeled  but  appear  in  three  distinct  groups  in  which 
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frequency  increases  from  left  to  right  in  steps  according  to  the  frequency  in¬ 
crement  chosen.  These  quantities  {gain,  phase)  are  calculated  only  if  SKIP  =  0. 
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IV.  The  Mechanics  of  Using  the  Program 


The  mechanics  of  using  the  program  are  very  simple.  They  consist  essen¬ 
tially  of  knowing  how  to  get  '''e  input  parameters  into  the  computer.  This  is 
accomplished  by  punching  the  values  of  the  input  parameters  on  a  set  of  cards 

5,'€ 

in  the  format  and  sequence  indicated  below,  placing  these  at  the  end  of  the 
program  deck,  and  instructing  the  operator  to  clear  the  core  and  load  the 
cards.  The  card  format  for  input  parameters  is  as  follows: 

Card  Format  for  Input  Parameters 

Card  1 

Column  1  requires  a  "1"  punch  if  it  is  desired  to  skip  to  a  new  page  for 
each  tube.  Columns  2-72  may  contain  any  identifying  information  for  the 
printed  output  such  as  tube  name  and  date. 

Card  2 


Parameter 

N  NCVS 

VO 

FIO  FO 

A 

B 

C 

Columns 

1-5  6-10 

11-18 

19-26  27-34 

35-41 

42-48 

49-55 

Type  Conversion 

I  I 

E 

E  E 

F 

F 

F 

Decimal  Places 

0 

0  0 

4 

4 

4 

Card  3 

Parameter 

NKK 

SFRIQ 

FFDEL 

Columns 

1  -5 

6-15 

16-25 

Type  Conversion 

I 

F 

F 

Decimal  Places 

4 

4 

=■'  See  FORTRAN  manual  for  an  explanation  of  conversion  types. 


Card  4 


Parameter 

DGAP(l) 

DGAP(2) 

DGAP(3) 

DGAP(4) 

DGAP(5) 

DGAP(6) 

Columns 

1-8 

9-16 

17-24 

25-32 

33-40 

41-48 

Type  Conversion 

F 

F 

F 

F 

F 

F 

Decimal  Places 

4 

4 

4 

4 

4 

4 

Card  5 

Parameter 

RQ(1) 

RQ(2) 

RQ(3) 

RQ(4) 

RQ(5) 

RQ(6) 

Columns 

1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

Type  Conversion 

F 

F 

F 

F 

F 

F 

Decimal  Places 

2 

•  2 

2 

2 

2 

2 

Card  6 

Parameter 

QO(l) 

QO{2) 

QO(3) 

QO(4) 

QO(5) 

QO(6) 

Columns 

1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

Type  Conversion 

F 

F 

F 

F 

F 

F 

Decimal  Places 

2 

2 

2 

2 

2 

2 

Card  7 

Parameter 

QEXT(l) 

QEXT(2) 

QEXT(3) 

QEXT(4) 

QEXT(5) 

QEXT(6) 

Columns 

1-10 

11-20 

21-30 

31-40 

41-50 

51-60 

Type  Conversion 

F 

F 

F 

F 

F 

F 

Decimal  Places 

2 

2 

2 

2 

2 

2 

Card  8 

Parameter 

S(l) 

S(2) 

S(3) 

S(4) 

S(5) 

Columns 

1-8 

9-16 

17-24 

25-32 

33-40 

Type  Conversion 

F 

F 

F 

F 

F 

Decimal  Places 

4 

4 

4 

4 

4 
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Card  9 


Parameter 

THETA 

Columns 

1 

-  10 

Type  Conversion 

F 

Decimal  Places 

4 

Card  10 

Parameter 

POTD 

WOT 

SKIP 

Columns 

1  -  10 

11  -  20 

21  -  30 

Type  Conversion 

F 

F 

F 

Decimal  Places 

2 

2 

2 

The  input  should  end  with  card  10  if  tunings  are  not  to  be  included.  If 
there  is  another  tube  to  be  analyzed,  the  next  card  will  be  of  the  same  format 
as  Card  1  above,  etc.  The  program  stops  when  it  runs  out  of  cards.  In  this 
case,  the  value  of  SKIP  must  be  nonzero.  Cards  11  and  12  are  to  be  included 
if,  and  only  if,  SI^P  =  0. 

Card  11 

Parameter  FNU(l)  FNU(2)  FNU(3)  FNU(4)  FNU(5)  FNU(6) 

Columns  1  -  10  11  -  20  21  -  30  31  -  40  41  -  50  51  -  60 

Type  Conversion  F  F  F  F  F  F 

Decimal  Places  5  5  5  5  5  5 

This  card  may  be  followed  by  as  many  cards  of  the  same  format  as  there 
are  sets  of  tunings  to  be  investigated. 

Card  12 

Punching  999-  99999''  in  columns  1-10  indicates  to  the  program  that 
there  are  no  more  sets  of  tunings  to  be  investigated  for  this  particular  tube 
and  to  proceed  to  the  next  tube. 
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MASTER  CONTROL  PROGRAM 


C 

C 


101 


107 


no 


REVISED  MCK  PROGRAM 
PART  1 

DIMENSION  DBGAINI 200) .REL( 10) 

DIMENSION  DGAP(IO) .RQ ( 10 ) »00{ 1 0 ) .OEXT ( 10) .S ( 10 ) . VI ( 10 ) »VOC( 10) .VS( 
110)  *RHO(  10)  .BP  ( 10)  .VBAR(  10)  .BE  ( 10  )  .FK(  10)  .BOPOSI  10.10)  .WQWPPdO)  .W 
20WP(10) .GPOS( 10) .BPOS(IO) .TPOS ( 10 ) . ZET AP ( 10. 10 ) .YOGOPI 10.10) .BMINU 
3S( 10) .GMINUSI 10 ) . TMINUS ( 10 ) .BOMIN ( 10. 10) .WOWPM( 1 0) .WOWM ( 10 ) . ZET AM ( 
A 10. 10) .YOGOM( 10. 10).EPOS( 1C ) .FPOS( 10) .CKP( 10.10) . FMP ( 10 . 10 ) . FHP ( 10 
5.1C) .EMIN( 10) .FMIN( 10).CKM( 10.10) . FMM ( 1 0. 1 0 ) . FHM ( 10. 1 0 ) .GELGOI 10) . 
6GELG0P( 1C ) .OEL( 10) .G( 10) .G0( 10) .GP( 10) .GM( 10) .D( 10) .QOEL(IO) .GPST( 
710 .GMST ( 1C ) .GAIN! 2C0) .PHASE! 20C ) .VPTR( 10) .VSUBD( 10) .FFO(200) 
DIMENSION  FNU( 10) .FREAL ( 10 . 1 0 ) .F 1  MAG ( 10.10) .Y ( 10 ) .GEL( 10) .ZGR( 10) . 
IZGI ( 10) .PPLUSR(20.20) .PPLUSI (20.20) .COLMR( 20.20) .COLMI (20.20) .FLMP 
2R(20.20 ) .FLMPI (20.20) .DR(2C.2C) .DI ( 20.20) .FIR(20.20) .FI  I (20.20) .FM 
3ULTR( 20.2C ) .FMULTI (20.20)  .SUMR ( 20 .20  )  .SUM I ( 20 . 20 ) .CR ( 20 . 20 ) . C I ( 20 . 
420) .ZGGLR(20.20) . ZGGL I ( 20 .20 ) 

DIMENSION  BEF( 3 . 1 C  .  1 0 ) . FKKKK ( 3 . 10 . 1 0 ) .GPOSS ( 3 . 1 0. 10) .BPOSS ( 3 ♦ 10. 10 
1 ) .TPOSS( 3.10. 10) .BQPOSS(3.10.10) .WQWPPP ( 3 . 10. 10 ) .WQWPW( 3. 10. 10 ) .ZE 
2TAPP( 3* 10.10) .Y0G0PP( 3.10.10) .GM INSS ( 3 . 10. 10 ) .BMINSS( 3. 10» 10) .TMlN 
3SS( 3.10.10) .B0MINN(3.10.10) .WQWPMM( 3 . 1 0. 10 ) .WCWMW( 3.10.10) .ZETAMM( 

43.10.10)  .YOGOMM( 3.10.10) .GELGPP ( 3 . 10 . 10 ) . EPOSS ( 3 . 10. 10 ) .FPOSS(3.10 

5.10)  .C<PP( 3.10.10) .FMPP(3.10.10) .FHPP( 3. 10. 10) .EM1NN( 3. 10. 10) .FMIN 
6N(3.10. 10) .CKMM( 3 .10.10) .FMMM( 3. 10. 10) .FHMM( 3 . 10. 10) .GELGOO( 3 ♦ 10. 1 
70) 

COMMON  BEE. FKKKK. GPOSS. BPOSS. TPOSS.BQPOSS. WQWPPP. WQWPW. ZET APP.YOGO 
IPP.GMINSS.BMINSS. TM INSS. BOM  INN. WOWPMM.WQWMW.ZE TAMM. YOGOMM. GELGPP. E 

2poss.fposs.ckpp.fmpp.fhpp.eminn.fminn.ckmm.fmmm,fhmm.gelgoo 


READ  INPUT  TAPE 
READ  INPUT  TAPE 
READ  INPUT  TAPE 
READ  INPUT  TAPE 
READ  INPUT  TAPE 
READ  INPUT  TAPE 
READ  INPUT  TAPE 
frio=sfrio 
NK  =  3 


2.1 

2 .2.  N. NCVS.VO.fi O.FO. A. B.C 
2 .9.NKI<^.SFRI0.FFDEL 

2.3.  (DGAP( I ) . 1=1 .NCVS) 

2.4.  (R0( I ) . 1=1 .NCVS  ) 

2.4.  (Q0( I ) .1=1 .NCVS  ) 

2 .4.  (OEXT( i ) . 1=1 .NCVS) 


FREOM=-0,0P00 

FDEL=O,O0OO 

NCVSL1=NCVS-1 

NCVSP1=NCVS+1 

READ  INPUT  tape  2 . 3 . ( S( I)  .  I  =  1  .NCVSL 1  ) 

READ  INPUT  TAPE  2. 5. THETA 

READ  INPUT  TAPE  2 .20 .POTD .WOT .SK I P 

EM=1,759E11 

CL=2.998E0 

EMC2=EM/CL**2 

THETA=THETA*3. 14159265/100.0 


FPSO=R,R54E-l? 

PI=3. 14159265 
IF  (POTD)  106.107.106 

CALL  POTDEP(NCVS.VO.FIO.A.B.C.DGAP.VPTR.RSUBD.VBAR.VSUBD) 
DO  no  1=1. NCVS 
VI ( I ) =VBAR( I ) 

G0( I ) =F lO/VI ( I ) 

CONTINUE 

VI (NCVSPl )=VSUBO(NrVSPl ) 


N.B.:  The  Bessel  function  subroutine  called  for  in  this  program  as  BESSELF  is  GEBSL,  Share 
distribution  number  271- 
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G0(NCVSP1 )=FI0/V1 (NCVSPl ) 

GO  TO  111 

106  DO  104  I=1,NCVSP] 

VI  (  I )  =VO 

104  GO( I)=FI0/VI(  I  ) 

111  WRITE  OUTPUT  TAPE  3.1 

WRITE  OUTPUT  TAPE  3 .6 .NCVS.N , VO 
WRITE  OUTPUT  TAPE  3 . 14.F 10. FO , A, B ,C 
WRITE  OUTPUT  TAPE  3.13 

WRITE  OUTPUT  TAPE  3 . 12 . ( V 1 ( I ) . I = 1 .NCVSPl ) 

WRITE  OUTPUT  TAPE  3 . 1 5 . ( DGAP ( I  ) .  l  =1  .NCVS ) 

WRITE  OUTPUT  TAPE  3. 16. ( R0( I ) . I = 1 .NCVS ) 

WRITE  OUTPUT  TAPE  3 . 1 7 . ( 00 ( 1) . I = 1 .NCVS ) 

WRITE  OUTPUT  TAPE  3. 1 8 . ( OEXT ( I  ) .  I  =  1 .NCVS ) 

WRITE  OUTPUT  TAPE  3. 19 . ( S( 1) . I =1 .NCVSLl ) 

WRITE  OUTPUT  TAPE  3 . 56 .NKK.FR 1 0. FFDEL 
IF(NKK-l)  320.321.322 

320  PAUSE 

321  NK=1.00 
FREQMrO.OO 

322  DO  103  K=1 .NK 
FFO(K )=FPEOM*FO+FO 
DO  115  1=1. NCVSPl 

VOC ( I ) =SORTF ( 1 .0-1 ,^/ ( 1 .0+EMC2*V I ( I ) ) ) 

VS( I ) =VOC( I )*CL 

RHO( I )  =  ( 39.37**2 )*FI0/(PI*(B**2-C**2)*VS(  I ) ) 

BP( I )=SQRTF (EM*RHO( I ) /EPSO*SQRTF ( (1.0-V0C( I )**2)**3) )/(VS( I )*39.37 

BE( I ) -2.0*PI*FF0(K)/(VS( I )*39,37 ) 

FK ( I )=BE ( I ) *VOC( I ) 

IF  (K-1)  120.121.122 

120  PAUSE 

121  GPOS( I ) =0.6*BE(  I  ) 

GMINUSI I ) =l.fl*RE(  I  ) 

GO  TO  115 

122  GPOS( I ) =0.97*GPST ( I ) 

GMINUSI I  ) =2.30*GMST( I  ) 

115  CONTINUE 
128  DO  105  L=1 .N 

DO  102  1=1. NCVSPl 
IF  (I-l)  160.161.162 

160  PAUSE 
162  111=1-1 

DO  163  I  1  =  1 .  1 1  I 
IF  (VI ( I )-VI ( I  I )  )  163.164.163 
163  CONTINUE 

161  GP(n=GPOS(I) 

GM( I )=GMINUS(  I  ) 

153  DDG=0.001 

154  CALL  SPCRDPI I .A.B.C.DDG.GP.BP.BE.FK.GGPOS.BBPOS.TTPOS) 

GPOS( I ) =GGPOS 

BPOS( I ) =BBPOS 
TPOS( I ) =TTPOS 
BOPOSIL. I )=BE( I )-BPOS( I ) 

WOWPP( I ) =BOPOS(L. I )/BP( I ) 

WOWP( I ) =BOPOS(L.I ) /BE(  I  ) 

ZETAPIL. I ) =BPOS( I )*BQPOS( L.I )*TPOS( I ) **2 / ( GPOS ( I )**2*(GP0S( I )**2+T 
1P0S( I )**2 ) ) 
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YOGOP (L, I )=BE  ( I ) / (BQPOS(L»l ) • ( 1 . 0  +  EMC2*V K I ) ) * ( 2 .0+EMC2 *VI  (  I  ) ) ) 
CALL  SPCROM( I , A , B * C * DOG *GM »BP ♦ BE ♦ FK » GGM I N* BBM I N * T TMl N ) 

GM1NUS( 1 )=GGMIN 
BM1NUS( 1  )=BBM1N 
TM INUS ( 1 )=TTM1N 
BQM1N(L,  n=-(BE(  I  )-BMINUS(  1  )  ) 

WQWPM( I ) =BQMIN( L* I ) /BP ( 1 ) 

WQWM( 1 )=PQMIN(L» I ) /BE(  I  ) 

ZETAM(L,  I  )=BMINUS( I  )*BQMIN(L.I )»TMINUS( I ) **2 / ( GM I NUS ( 1) **2*  (  GM 1 NUS 
1 ( 1 )**2+TMINUS( 1 )**?) ) 

YOGOM(L»  U  =BE( I ) / (BQMIN(L» I  »»(! .0  +  EMC2*VI ( I ) ) * ( 2 .0  +  EMC2*Vl ( 1) ) ) 

GO  TO  181 

164  GPOS( 1 ) =GPOS( 11 ) 

BPOS( I )=BPOS(  I  I  ) 

TPOS( I ) =TPOS( I  I ) 

BOPOS(L»I )=BOPOS(L,I  I  ) 

WQWPP( 1 ) =WQWPP(  I  I  ) 

WQWP( 1 1 =WQWP(  I  I  ) 

ZETAP  (L,  I  )  =Z'"TAP(  L.I  I  ) 

YOGOP(L, I )=YOGOP(L»I  I  ) 

GMINUS( I  )=GMINUS(  I  I  ) 

PM1NUS( I  )=BMINUS(  I  I  ) 

TMINUSC I )=TMINUS(  I  I  ) 

BOMIN(L,  I)=BOMIN(L,n  J 
WQWPM( I ) =WOWPM( I  I  ) 

WQWM(  I ) =WQWM(  I  I  ) 

ZETAM(L, I  )  =ZETAM( L ,11) 

YOGOM(L  ,  I  )  =YOGOM.(  L,  I  I  ) 

181  BEE(K,L,  n=BF(  I  ) 

FKKKK  (K,L,  I  )='='<(  I  ) 

GPOSS(<,L, I )=GPOS( I ) 

OPOSS(K,L, I )=BPCS(  I  ) 

TPOSS(<,L, I )=TPOSJ  I  ) 

BQPOSS(K,L,I )=BQPOS(L,I ) 

WQWPPP( K ,L , I ) rWQWPP ( t ) 

WQWPW(K,L, I )=WQWP ( I ) 

ZETAPP(K,L,I  )=ZFTAP(L,I ) 

YnGOPP(K.L.I )=YOGOD(L,I  ) 

GMINSS(K,L.I )=GMINUS(  I  ) 

PMINSS(K,L.I  )=BVIINUS«  I  ) 

TMIN5S(K,L, I ) =TMINUF  (  I  ) 

BQMINN(K,L,I )=BOMIN(L,I  ) 

WQWPMM( K ,L , I ) =WQWPM(  I  ) 

WQWMW ( K , L , I ) =WOWM ( I ) 

ZETAMM( K .L , I)=ZFTAM(L,I) 

YOGOMM(K,L,I )=YOGOM(L,I  ) 

IF  (L-1)  12E,126»10? 

125  PAUSE 

126  GPST( I )=GDnS( I ) 

GMST( I ) =GMINUS(  I  ) 

102  CONTINUE 

CALL  EFCMH(L,A,d*C,NCVS,V  I  ,TPOS,GPOS»BPOS»EPOS»FPOS,CI^P  ,DUM1  ,DUM2* 
10UM3) 

CALL  EFCMH(L,A,B,C,NCVS,VI,TMINUS,GMINUS»BMINUS,EMIN,FMIN,C<M,DUM1 
1 ,DUM2 ,DUM3) 

DO  108  I=1,NCVS 

FMP( L  ,  I  ) =BESSELF ( DGAP( I ) *BPOS( I )  ,0, 1 ) 

FMM( L  ,  I  )  =BESSELF ( DGAP ( I ) *BMINUS(  I  )  ,0, 1 ) 


16 


108 

105 

103 


175 

400 

176 


189 

179 
401 

180 

323 

325 


FHP(L.  I  )=EPOS(  I  )  +  (ZETAP(L»n*FPOS(  I  )  ) 

FHM(L. I )=EMIN( I ) - ( ZETAM( L* 1) *FM IN ( I ) ) 

GELG0(  I  )=0.25*(FHP(L»  I  )»Y0G0P  ( L»  I )  *CKP  (  L  ♦  I  )  ♦♦2*FMP  (  L  ♦  I  )♦^^2/(  ( 1.0+Z 
1ETAP(L» I ) )»»2)“(FHM(L* I )*Y0G0M(Lt I ) *CKM ( L ♦ I ) ♦♦Z^FMM ( L ♦ I ) **2/ ( (1.0- 
2ZETAM(L»I ) )**2) ) } 

GELG00(K,L. I )=GELG0( I ) 

EP0SS(K»L. I )*EP0S( I ) 

FP0SS(K»L.I )=FP0S(  I  ) 

CKPP(K.L.I  )«C<P(L.n 
FMPP(K»L» I )»FMP(L«I  ) 

FHPP(K.L. I )«FHP(L»I  ) 

EMINN(K»L. I )*EVIN( I ) 

FMINN(K»L. I )*FMIN( I ) 

CKMM(K.L» I )«CKM(L» I ) 

FMMM(K»L»I )«FMM(L«I) 

FHMM(K.L» I )*FHM(L«I  ) 

CONTINUE 

CONTINUE 

FREOMrFREQM+FDEL 

CONTINUE 

MM  =  0 

IF(SKIP)  176.175.176 

READ  INPUT  TAPE  2 .49 ♦ ( FNU( 1) » I »1 »NCVS ) 

IF(FNU( 1 )-999.999)  400»101»101 

FRIQ»SFRIQ 

DO  171  K»1«NKK 

FF0(K)=FRI0*F0+F0 

DO  189  Isl.NCVS 

GELGOP( I )=0»0 

CONTINUE 

IF  (WOT)  179.180.179 

IF(MM)  180.401.180 

WRITE  OUTPUT  TAPE  3. 10. FRIO 

DO  172  L=1.N 

IF(NKK-l)  323.324.325 

PAUSE 

CALL  TERPOL(NCVSPl.L.BEE.FRIO.BE) 

CALL  TERPOL (NCVSPl.L.BEE.FRIO.BE) 

CALL  TERPOL(NCVSPl.L.FKKKK.FRIO.FK) 

CALL  TERPOL ( NCVSPl .L .GPOSS .FR lO.GPOS ) 

CALL  TERPOL ( NCVSPl .L .BPOSS .FR 10 .BPOS ) 

CALL  TERPOL(NCVSPl.L.TPOSS.FRIQ.T?OS) 

CALL  TERPOL (NCVSPl .L .WQWPW.FR lO.WOWP ) 

CALL  TERPOL ( NCVSPl. L.WQWNW.FRIQ.WOWM) 

CALL  TERPOL ( NCVSPl .L .WQWPPP.FR lO.WOWPP ) 

CALL  TERPOL (NCVSPI.L.GMINSS.FRIQ.GM I NUS) 

CALL  TERPCl (NCVSPl .L. BMIMSS.fr I Q.BM I NUS ) 

CALL  TERPOL (NCVSPl.L .TMINSS.FRIO.TM INUS ) 

CALL  TERPOL (NCVSPl.L.WOWPMM.FRin.WQWPM) 

CALL  A  INTER (NCVSPl .L .BGPOSS.FR  lO.BOPOS  ) 

CALL  AINTER(NCVSP1.L .ZETAPP.FRIQ.ZETAP) 

CALL  A  INTER (NCVSPl.L . YOGOPP . FR 1 0 . YOGOP ) 

CALL  A  INTER (NCVSPl.L .BQM I NN . FR 1 0 . BQMI N ) 

CALL  AINTER (NCVSPl.L .ZETAMM.fr 10. ZET AM) 

CALL  AINTER ( NCVSPl.L. YOGOMM.fr 10. YOGOM) 

CALL  TERPOL (NC VS. L.GELGOO. FRIO. GELGO) 
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324 


349 


350 

351 

183 

178 

402 


CALL  TERPOL(NCVS.L.EPOSS»FRIQ*EPOS) 
CALL  TERPOL(NCVS»L.FPOSS»FRIO*FPOS) 
CALL  TERP0L(NCVS*L»EM1NN  »FR10*EMIN) 
CALL  TERP0L(NCVS»L*FM1NN  tFRIQfFMIN) 
CALL  A1NTER(NCVS»L»CKPP»FRI0»CKP) 
CALL  AINTER(NCVS.L*FMPP»FRIQ«FMP) 
CALL  A  INI ER(NCVS»L»FHPP*FR10»FHP) 
CALL  AINTER(NCVS.L*CKMM»FR10.CKM) 
CALL  AINTER(NCVS.L»FMMM»FRIO»FMM) 
CALL  A1NTER(NCVS»L«FHMM*FRI0»FHM) 

GO  TO  351 

DO  349  1=1*NCVSP1 

BE(  1  )=BEE( 1*L» 1 ) 

FK( I)=FK<K<(1«L«1 ) 

GPOS( I )=GP0SS(  1 «L«1  ) 

BPOS( 1 )=BP0SS( 1 »L»1 ) 

TPOS( I )=TP0SS( 1«L»1 ) 

WOWP(  I  ):WOWPW(  1  «L«n 
WGWM( I )=WOWMW(l»L«l 1 
WOWPP ( I ) sWOWPPP ( 1 «L  « 1) 

GM1NUS( 1 )=GM1NSS( liL*I ) 

BM1NUS( I )=BM1NSS( 1«L*I ) 

TMINUS( 1 )=TMINSS( 1*L#I ) 

WOWPM ( I ) *WOWPMM ( 1  *  L  » 1) 

BOPOS(L« 1 )=BOPOSS( 1*L» I ) 

ZETAP(L«1 )»ZETAPP(1«L»I) 

YOGOP(L» 1 )=YOGOPP( 1*L» 1 ) 

B0M1N(L* 1 )*B0M1NN( 1«L» I ) 
ZETAM(L»1)=ZETAMM(1,L»I) 

YOGOM(L» 1 )aYOGOMM( 1 ^L* I ) 

CONTINUE 
DO  350  1*1»NCVS 
GELGO(  1  )=GELGOO(l»L»n 
EPOS(l)*EPOSS(l»L»l) 

FPOS( I )«FPOSS( 1 .Ltl) 

EM1N( I )=EM1NN( 1 »Lil ) 

FM1N( I )aFMINN( 1 »L»I ) 

CKP(L» 1 )=CKPP( 1 tL*! ) 
FMP(L»1)*FMPP(1»L»1) 

FHP(L» 1 )=FhPP( 1 tLtl ) 

CKM(L»1 )=CKMM( 1»L»1  ) 

FMM(Lf  1  )=FMMyi(  1  tL*!  ) 

FHM(L» 1 )=FHMM( 1 tL*!  ) 

CONTINUE 

DO  183  I=1»NCVS 

GELGOP( I )=GELG0P( I )+GELGO( I ) 


CONTINUE 

IF  (WOT)  178»172»178 
IF(MM)  172.402»17? 

WRITE  OUTPUT  TAPE  3fll»L 

WRITE  OUTPUT  TAPE  3»26 ♦ ( REI I ) » 1=1 »NCVSP1 ) 

WRITE  OUTPUT  TAPE  3 tT* ( BPOS( I ) » 1=1 iNCVSPl ) 
WRITE  OUTPUT  TAPE  3 *8 ♦ ( BMINUS ( I ) ♦ I = 1 iNCVSPl ) 
WRITE  OUTPUT  TAPE  3 » 2 1  ♦  (  GPOS  (  I  )  ♦  I  =  1  »NC\/SP1 ) 
WRITE  OUTPUT  TAPE  3 *22 » ( GMI NUS ( I ) » I =1  * NCVSPl ) 
WRITE  OUTPUT  TAPE  3 *23  * ( TPOS (I ) *  I  =  1  * NCVSP 1 ) 
WRITE  OUTPUT  TAPE  3 *24 ♦ ( TMINUS ( I ) » I =1 »NCVSP1 ) 
WRITE  OUTPUT  TAPE  3 * 66 . ( BQPOS ( L *  I ) » I  =  1 ♦ NCVSP 1 ) 
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WRITE  OUTPUT  TAPE  3  *  67 » ( BOMIN ( L * L) *  I  =  1  * NCV5P 1 ) 

WRITE  OUTPUT  TAPE  3 , 25 » ( WOWPP ( I ) » I = 1 tNCVSPl ) 

WRITE  OUTPUT  TAPE  3 ♦ 27  ♦  ( WOWPM ( I ) *  I =1 tNCVSP 1 ) 

WRITE  OUTPUT  TAPE  3 ♦ 28 ♦ ( WOWP ( I ) *  I  =  1 tNCVSPl ) 

WRITE  OUTPUT  TAPE  3 *  29  ♦  ( WQWM( I ) *  I  =  1 tNCVSPl ) 

WRITE  OUTPUT  TAPE  3 ♦ 30  ♦  ( ZETAP ( L *  I ) *  I  =  1 *NCVSP 1 ) 

WRITE  OUTPUT  TAPE  3 ♦ 31  * ( ZETAM ( L » I ) » I « 1  * NCVSP 1 ) 

WRITE  OUTPUT  TAPE  3» 32 ♦ ( YOGOP < L *  I ) ♦ I  =  1 »NCVSP 1 ) 

WRITE  OUTPUT  TAPE  3 ♦ 33  * ( YOGOM ( L » I ) *  I  =  1 »NCV5P 1 ) 

WRITE  OUTPUT  TAPE  3* 34 ♦ ( EPOS ( I ) *  I  =  1 •NCVS ) 

WRITE  OUTPUT  TAPE  3 » 35  ♦  ( EMI N ( I ) *  I  =  1  * NCVS ) 

WRITE  OUTPUT  TAPE  3 ♦ 36 ♦ ( EPOS ( I ) *  1  =  1 »MCVS ) 

WRITE  OUTPUT  TAPE  3» 37 ♦ ( FMIN ( 1) ♦ I *1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 ♦ 38 ♦ ( CKP ( L ♦ I ) ♦ I  * 1 ♦ NCVS ) 

WRITE  OUTPUT  TAPE  3 ♦ 39  ♦  ( CKM ( L » I ) ♦ I =1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 tAO ♦ ( FMP ( L ♦ I ) ♦ I =1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3^41 ♦ ( FMM(L ♦ I ) ♦ I =1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 ♦ 42 ♦ ( FHP ( L ♦ I ) ♦ I =1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 ^43 ♦ ( FHM( L ♦ 1) ♦ I =1 ♦ NCVS ) 

WRITE  OUTPUT  TAPE  3 ♦ 44 ♦ ( GELGO ( I ) ♦ I =1 ♦NCVS ) 

172  CONTINUE 

IF  (WOT)  184,185tl84 

184  IF(MM)  185^403tl85 

403  WRITE  OUTPUT  TAPE  3t 45 ♦ ( GFLGOP ( I ) ♦ I »1 ♦NCVS ) 

DO  118  l=ltNCVS 

OEL(  n=1.0/(PO(  1  )*GELGOP(n*GO(  I  )  ) 

0( I ) =1 #0/ ( 1 #0/00 ( I )+l«0/OEL( 1 ) +1»0/OEXT ( I ) ) 

OOEL(  n  =  l,0/(1.0/0(I  )-l  fO/OEXT  (  I  )  ) 

GEL (I )=GELGOP( 1 )*G0( I ) 

REL( I )=1.0/GEL( I ) 

118  CONTINUE 

WRITE  OUTPUT  TAPE  3 ♦AS  ♦  ( QEL ( I  )  ♦  I  =  1 tNCVS ) 

WRITE  OUTPUT  TAPE  3^47  ♦  (0( I ) ♦ I *1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 ^48  ♦  ( OOEL ( I ) ♦ I *1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 ♦ 53  ♦  ( GEL ( I )  ♦  I  =  1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3  ♦  81  ♦  ( REL ( I  )  ♦  I  =  1 ♦ NCVS ) 

185  IF(SKIP)  177^186^177 

186  NN=2*N 
DELO=(FFO(K)-FO)/FO 

CALL  FREQIO^FNU^DELO^NCVS^FREAL^FIMAG ) 

DO  143  1=1  ♦NCVS 
Y( I  )=Q( I )/OEL( I ) 

ZGR(I)=Y(I)*FREAL(I  )/GEL(I) 

ZGI(I)=Y(I)*FIMAG(I)/GEL(I) 

CALL  FOMATR ( N^ I ♦NCVS ♦CKP  tCYM ♦ FMP ♦ FMM ♦ FHP ♦ FHM ♦ ZETAP ♦ZET AM ♦ YOGOP ♦YOG 
lOM^PPLUSR^PPLUSI ♦COLMRtCOLfl ♦FLMPR»FLMP I ) 

140  CALL  DMATRX ( NCVSPl ♦NN *  I *BOPOS ♦BQMI N ♦BL ♦S ♦DR ♦D I ♦F I R ♦F I  I  ) 

IF  (I-l)  130^131^132 

130  PAUSE 

132  IF  (I-NCVS)  133^134^130 

131  CALL  CMMP ( DR^DI ♦NN*NN*COLMR^COLMI ,1 ♦FMULTR^FMULTI  ) 

GO  TO  143 

133  DO  170  JJ=1^NN 
DO  170  KK=1^NN 

ZGGLR  (  JJ^KK)  =G0(  I  )^'ZGR  (I  )<<^FLMPR(  JJ,ICK)-GO(  I  )*ZGI  (  I  )*FLMPI  ( JJ^KK) 
ZGGLI ( JJ*KK )=G0( I )*ZGI (I ) *FLMPR ( J J *KK ) +G0 ( I )*ZGR( I )*FLMPI ( JJ*KK) 
170  CONTINUE 

CALL  CMSUBT(ZGGLR^FIR^NN^SUMR  ) 
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CALL  CMSUBT (ZGGLI *FI I .NN.SUMI )  . 

CALL  CMMP  (  S'JMR  tSUM  I  * NN  *NN  »FMULTR  tFMULT  I  *1  *CR  *CI  ) 
DO  141  JJ=1»NN 
KK=  1 

FMULTR(JJ*KK)=CR(JJ*KK) 

FMULTI (JJ*KK)=CI ( JJtKK  J 
141  CONTINUE  .  , 

CALL  CMMP (DR *01 ,NN »NN ♦ FMULTR » FMULT I » 1 *CR *C I ) 

DO  142  JJ=1»NN 


142 
134 

143 


155 

156 
158 

157 


173 

404 


177 

171 


KK=  1 

FMULTR(JJiKK)=CR( JJ*KK  ) 

FMULTI (JJ*KK)=CI ( JJ»KK  ) 

CONTINUE 

CALL*^CMMP  (PPLUSR«PPLUS  I « 1  »NN*FMULTR  ♦FMULT  I  *1  ♦CR*C1  ) 

FIGVGR=CR(lil) 

FIGVGI*CI (1*1) 

CONTINUE 

EXPR=COSF(THETA) 

EXPI=SINF(THETA) 

FIVR=FIGVGR*GO(NCVS) 

FIVI  =  FIGVGI*GO(NCVS)  _ _ _ 

F  IVERsFREAL  (NCVS  )*F  I\/R*EXPR-F  IMAG  ( NCVS )  *F  I  VR^EXPI-FREAL  (  NCVS  )  *FI  VI 

1 *Fypi -F I  MAG (NCVS)*F1VI*EXPR 

FIVEI=FREAL(NCVS)*FIVR*EXPI+FIMAG(NCVS)*FIVR*EXPR+FREAL(NCVS)*FIVI 

1*EXPR-FIMAG(NCVS)*FIVI*EXP1 

GVR=- (RQ(NCVS)*QEXT( NCVS ) *F IVER ) / ( 1 »0+OEXT (NCVS) /OOEL( NCVS ) ) 
GVI=-(RQ(NCVS )*QEXT (NCVS )*FIVEI )/( 1 .O+OEXT ( NCVS ) /OOEL ( NCVS ) ) 

PHASE (K)=ATANF(GVI /GVR ) 

IF  (GVR)  155*156*156 
PHASE (K)sPHASE(K)+PI 
GO  TO  157 

IF  (PHASE(K))  158*157*157 
PHASE(K)=PHASE(K)+2.0*PI 
FONE=FREAL( 1)**2+FIMAG( 1 )**2 
FTW0=FREAL(NCVS)**2+FIMAG(NCVS)**2 

FIGV=FIVR**2+FIVI**2  ^  ^ 

GAIN(K)=4.0*RQ( 1 )*RQ(NCVS)*0EXT(1 )«0EXT(NCVS)*F0NE*FTW0*FIGV/( (1.0 
1+QEXT ( 1 ) /OOEL ( 1 ) )**2*( l.O+OEXT (NCVS) /OOEL (NCVS) )**2 ) 
0BGAIN(K)=4.3429*L0GF;gAIN(IC)  ) 

IF  (WOT)  173*177*173 
IF(MM)  177*404.177 

WRITE  OUTPUT  TAPE  3 » 52  * ( Y ( I ) *  I  =  1 *NCVS ) 

WRITE  OUTPUT  TAPE  3 » 54 . ( ZGR ( I ) ♦ I = 1 ♦NCVS ) 

WRITE  OUTPUT  TAPE  3 . 55  ♦  ( ZG I ( I ) *  I  =  1 *NCVS ) 

WRITE  OUTPUT  TAPE  3 *70 .F I VR tF I V I *GVR *GV I 
FRI0=FRIQ+FF0EL 
CONTINUE 

IF  (S)CIP)  101.187*101 

3.50.(FNU(I  ). 1  =  1. NCVS) 

3.71 *( GAIN! KK) *KK=1 *NKK) 

3.72. (DeGAIN(KK) *KK=1 *NKK) 

3.72 . (PHASE (KK) .KK=1 *NKK) 


FORMAT  STATEMENTS 
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WRITE 

OUTPUT 

TAPE 

WRITE 

OUTPUT 

TAPE 

WRITE 

OUTPUT 

TAPE 

WRITE 

OUTPUT 

TAPE 

MM=1 
GO  TO 

175 

C 

1 

FORMAT  (54H 

2 

FORMAT  (215. 

3E8.0 

20 


•3 

FORMAT 

( 7FBt4  ) 

u 

FORMAT 

(7F10.2) 

5 

FORMAT 

(Flo. 4) 

6 

FORMAT 

(20H  NUMBER  OF  CAVI TIES=I 3/17H  NUMBER  OF  M0DES*I3/12H  DC  VO 

1LTAGE  =  E10.3  ) 

7 

FORMAT 

(25H 

R  + 

»7F12.6) 

8 

format 

(25H 

R- 

.7F12.6) 

9 

FORMAT 

( 15 .2F10 

.4  ) 

10 

FORMAT 

(  /////'6H 

FDEL=F10.4//) 

1 1 

FORMAT 

(///3H  M 

=  15/) 

12 

format 

(25H 

DC  CAVITY  VOLTAGE 

.7E12.3) 

13 

FORMAT 

(7H  CAVITY26X.1H1 11X.1H211X. 

1H31 1X.1H411X.1H511X.1H611X»1H7 

14 

1  ) 

format 

( 12H  DC 

CURRENT=F10.3/30H 

CENTER  FREQUENCY! CYCLES/SEC )=E1 2. 

14/3H  Ar 

F7.4/3H 

prF7.4/3H  C=F7.4////) 

15 

FORMAT 

(25H 

D 

»7F12.4) 

16 

format 

(25H 

R/0 

.7F12.2) 

17 

FORMAT 

(25H 

00 

.7F12.2) 

18 

FORMAT 

(25H 

OEXT 

.7F12.2) 

19 

FORMAT 

(29H  CENTER  GAP  TO  CENTER 

GAP  »7F12.4) 

20 

FORMAT 

( 3F10.2) 

21 

FORMAT 

(25H 

G+ 

»7F12.6) 

22 

FORMAT 

(2  5H 

G- 

.7F12.6) 

23 

FORMAT 

(  25H 

T  + 

»7F12.6) 

24 

format 

(25H 

T- 

»7F12.6) 

25 

FORMAT 

(25H 

WO+/WP 

»7F12.6) 

26 

FORMAT 

(25H 

PE 

»7F12.6) 

27 

FORMAT 

(25H 

WO-/WP 

.7F12.6) 

28 

format 

(25M 

WO+/W 

»7F12.6) 

29 

FORMAT 

(25H 

WO-/W 

»7F12.6) 

30 

FORMAT 

(25H 

ZETA+ 

»7F12.6) 

31 

format 

(25h 

zeta- 

»7F12.6) 

32 

format 

(25H 

YO+/GO 

»7F12.6) 

33 

format 

(25H 

YO-/GO 

.7F12.6) 

34 

FOR^'AT 

(25H 

E  + 

.7F12.6) 

35 

format 

(25H 

E- 

.7F12.6) 

36 

format 

(2'^H 

F  + 

.7F12.6) 

37 

POPMAT 

(25H 

F- 

*7F12.6) 

38 

format 

(25H 

C  + 

»7F12.6) 

39 

FORMAT 

(  25H 

c- 

»7F12.6) 

40 

FORMA'" 

(25H 

M+ 

.7F12.6) 

41 

FORMAT 

(25H 

M- 

»7F12.6) 

42 

FORMAT 

(25H 

H  + 

»7F12.6) 

43 

format 

(25H 

H- 

»7F12.6) 

44 

format 

(25H 

GELGO 

.7F12.6) 

45 

format 

(/////25M  GELGOP 

*7F12.6) 

46 

FORMAT 

(25H 

OEL 

.7F12.6) 

47 

format 

(25H 

0 

.7F12.6) 

48 

format 

(25H 

QOEL 

.7F12.6) 

49 

format 

(7F10.5) 

50 

format 

(  //25H 

NU 

I7F12.6 ) 

52 

format 

(  25H 

Y 

I7F12.6) 

53 

format 

(25H 

GEL 

.7F12.6) 

54 

format 

(25H 

ZGR 

I7F12.6) 

55 

format 

(25H 

ZGT 

.7F12.6////) 

56 

format 

(5H  NKK 

--I5/6H  FREQ  =  F10.4/5H 

DEL=F10.4/ ) 

66 

format 

(  25H 

BQPOS 

.7F12.6) 

67 

format 

(  25H 

BOMIN 

.7F12.6) 

Z1 


70  FORMAT  (6H  FIVR  =  F10i6/6H  FI VI =F10.6/5H  GVR  =  F10»6/5H  GVI=F10»6) 

71  FORMAT  (//8ElA,6/( 8E14.6 ) ) 

72  FORMAT  ( //8 FI A ,6/ ( 8F 14 ,6  I ) 

81  FORMAT(25H  REl.  7F12.6) 

END 


22 


POTENTIAL  DEPRESSION 


SUBROUT INE  POTDEP (NCVS » VO *F I 0» A . B »C f DGAP» VPTR , RSUBD* VBAR .VSUBD ) 
DIMENSION  VOC( 10) .DGAPI 10) ♦VPTRI 10) *VSUBD( 10) .VBAK ( 10 ) *R0EPS( 10) .0 
IDDEVI 25  ) 

NCVSP1=NCVS+1 

EM=1.759E11 

CL=2.998E8 

EMC2=EM/CL**2 

VOCINCVSPl )=SQRTF( 1. 0-1. 0/(1. 0+EMC2*VO)*«2) 

ODDEVd  )=2.0 

0DDEV(2 )=2.0 

DO  A2  J=A.2A,2 

ODDEV ( J-1 ) =ODDEV ( J-2 ) +2 . 0 

ODDEVI J)=ODOEV( J-]  ) 

A?  CONTINUE 

Pl=3. 1A159265 
FMG0=376.7 
BPC=B**2+C**2 
BMC=B**2-C«»2 

VPTRINCVSPl )=1.0/ (2.0»PI ) *FMGO*F I O/VOC ( NCVSPl ) * ( 0. 5+LOGF ( A/0 ) -BPC/ 
1(A.O*BMC)+C««4*LOGF(B/C)/BMC»*2 ) 

VSUBD(NCVSPl)=VO-VPTR(NCVSPl)/( 1 . 0-VPTR ( NCVSP 1 )/(2.0*VO) ) 

1=0 
RD  =  B 

RRD=(B-C)/2.0 
GO  TO  31 

17  RD=RD-RRD 

31  1=1+1 

TEST=BMC/2.0-BPC/A.O+C»*A»LOGF(B/C)/BMC-( B**2-RD**2 )/2.0+C**2*LOGF 
KB/RD) 

IF  (I-l)  10.11.12 

10  PAUSE 

12  IF  (TEST)  13.14.15 

13  SIGN2=-1.0 
GO  TO  16 

15  SIGN2=+1.0 

16  IF  (SIGN+SIGN2)  17,18.17 

11  IF  (TEST)  19,14,20 

19  SIGN=-1.0 
GO  TO  17 

20  SIGN=+1.0 
GO  TO  17 

18  RDD=1.0E-6 

IF  (RRD-RDD)  14,14,21 

21  RD=RD+RRD 
RRD=RRD/10.0 
GO  TO  17 

14  RSUBD=RD 

DO  22  I=1,NCVS 

VOC( 1 )=SQRTF( 1.0-1.0/(1.0+EMC2*VO)«*2) 

ROEPS( I ) =FMGO/PI»FIO/VOC( 1 )»1,0/ ( BMC ) * ( 39. 37 ) **2 / ( 1 . 0-VPTR ( 1 )/ ( 2.0 
1*V0) ) 

M=0 

27  M=M+1 

FM  =  M 

IF  (FM/ODDEV(M)-1.0)  40,27,41 
41  PAUSE 

40  FPD  =  FM*P1/(2.0*DGAP(  1  )  ) 

IF  (C)  28,29,30 
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28  PAUSF 

30  FIlCD=BESSELF(FPD*Ctl.2) 

GO  TO  32 

29  FI1CD=0.0 

32  FK1BD=BESSELF( FPD*B.l .3) 

FI0RD=BESSELF(FPD*RSUBD.O.2( 

FK0RD*BESSELF(FPD*RSUBD»0.3 ) 

UM=( 1 .0/FPD)**2»( 1.0-FPD*B*FKlBD*FI0RD-FPD*C*FK0RD*FIlCD)/(39.37)# 
1*2 

ANR=8.0/(FM**2*PI**2»*ROEPSn ) *S INF ( FM*P I /2 .0 ) **2*UM/ ( 1 .0-ROEPS ( I ) 
1*UM/(2.0*V0) ) 

IF  (M-1 )  23.2At25 

23  PAUSE 

25  IF  (ANR-CHECK)  26.38»38 

24  CHE=1.0E-4 
CHECK=ANR*CHE 
SUM=ANR 

GO  TO  27 
38  SUM=SUM+ANR 
GO  TO  27 

26  VRAR(  n=VSUBD(NCVSPl  )-SUM 
2?  CONTINUE 

RETURN 

END 
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C  SPACE  CHARGE  REDUCTION  FACTOR 

SUBROUTINE  SPCRDP ( I . A tB *C .DDG *GP *  BP  *  BE . FK *GG* BB . T ) 

DIMENSION  GP( 10) .BP( 10) »BE( 10) »FK(10) 

G  =  GP(  I  ) 

J  =  0 

DG=DDG 
20  G  =  G+DG 

TTT=G**2»( BP( I )**2/( ( S0RTF(G»»2+FK( I )**2)-BE( I ) )»*2)-1.0) 

IF  (TTT)  20*20.12 
12  TT=SORTF(TTT) 

J  =  J+1 

IF  (C)  1.3.2 

1  PAUSE  11111 

2  FJOTC  =  BESSELF( TT*C.0.1  ) 

FJ1TC  =  BESSELF( TT*C.l .  1 ) 

FI0GC=BESSELF(G*C.0.2) 

FI 1GC=BESSELF(G*C.1.2) 

FNOTB=BESSELF( TT*B.n.A) 

FN1TB=BESSELF(TT*B.1.A) 

FNOTC=BESSELF( TT*C.O.A) 

FN1TC=BESSELF( TT*C.1.A) 

3  FI0GA=BESSELF(G*A.0.2) 

FI0GB=BESSELF(G*B.0.2 ) 

FI 1GB=BESSELF(G*B.1.2) 

FK0GA=BESSELF(G*A.0.3) 

FK0GB=BESSELF(G*B.0.3) 

FK1GB=BESSELF(G*B.1.3) 

FJOTB=BESSELF( TT^B.0.1 ) 

FJ1TB=BESSELF( TT*B.l  .1 ) 

BCOTH=(FKOGA*FI 1GB+FK1GB*F10GA )/(FKOGB*FIOGA-FKOGA*FIOGB) 

IF  to  A. 5. 6 
A  PAUSE  AAAAA 

5  FOFI=TT*B*FJ1TB-G*B*BCOTH*FJOTB 
GO  TO  7 

6  FOFI=( TT*B) **2*IFJ1TB*FN1TC-FJ1TC*FN1TB)-TT*B*G*B*BCOTH*(FJOTB*FN1 
1TC-FJ1TC*FN0TB)  +  TT»G*B**2^‘FI  1GC/FI0GC*(  FJ1TB*FN0TC-FJ0TC*FN1TB  )-(G 
2*B)**2*BC0TH*FI 1GC/FI0GC*(FJ0TB*FN0TC-FJ0TC*FN0TB) 

7  IF  (J-1)  A. 70. 69 

69  IF  (FOFI)  80.10.82 

80  SIGN2=-1.0 

GO  TO  83 

82  SIGN2=+1.0 

83  IF  (SIGN+SIGN2)  20.8.20 

70  IF  (FOFI)  6A. 10.66 
6A  SIGN=-1.0 

GO  TO  20 
66  SIGN=+1,0 

GO  TO  20 

8  nGG=1.0E-6 

IF  (DG-DGG)  10.10.9 

9  G=G-DG 
DG=DG/10,0 
GO  TO  20 

10  T=TT 
GG  =  G 

BB  =  SORTF (GG**2  +  FK( I )**2  ) 

RETURN 

END 
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c  SPACE  CHARGE  REDUCTION  FACTOR 

SUBROUTINE  SPCRDMI I , A *B * C  tDDG *GM *  BP ♦ BE  » FK » GG* BB . T ) 

DIMENSION  GM(  10)  !>BP(10)  .HE(  10)  .FK(IO) 

G  =  GM(  I  ) 

J  =  0 

DG=DDG 
20  G  =  G-DG 

TTT=G»»2*(BP( I )»*2/( (SORTF(G**2+FK( I ) **2 ) -BE ( I ) ) »»2 ) -1 . 0 ) 

IF  (TTT)  20.20,12 
12  TT=SQRTF (TTT) 

J  =  J  +  1 

IF  (C)  1,3,2 

1  PAUSE  11111 

2  FJ0TC=BESSELF(TT»C,0,1 ) 

FJ1TC=BESSELF(TT*C,1,1) 

FIOGC=BESSELF(G*C ,0,2 ) 

FI1GC=BESSELF(G*C,1,2) 

FNOTB=BESSELF( TT*B,0,A ) 

FN1TB=BESSELF(TT*B,1,A) 

FNOTC=BESSELF( TT*C,0,A) 

FN1TC=BESSELF(TT*C,1,A) 

3  FI0GA=BESSELF(G*A,0,2) 

FI0GB=BESSELF(G*B,0,2) 

FI1GB=BESSELF(G*B,1,2) 

FK0GA*BESSELF(G*A,0,3) 

FK0GB=BESSELF(G*B,0,3) 

FK1GB=BESSELF(G*B,1,3) 

FJ0TB=BESSELF(TT*B,0,1 ) 

FJ1TB=BESSELF(TT*B,1,1 ) 

BCOTH*  (  FK0GA*F  I IGB+FK  1GB*F  I  OGA  )  /  (  FK0GB*F  lOGA-FICOGA^F  1 0GB ) 

IF  (C)  4,5,6 
A  PAUSE  444A4 

5  fofi=tt*b*fjitb-g*b*bcoth*fjotb 

GO  TO  7 

6  FOFI=(TT*B)**2*(FJ1TB*FN1TC-FJ1TC*FN1TB)-TT*B»G*B*BCOTH*(FJOTB*FN1 
1TC-FJ1TC*FN0TB)+TT*G*B**2»FI 1GC/FI0GC*( FJ1TB*FNoTC-FJOTC*FN1TB )-(G 
2»B)»*2*BCOTH*FI1GC/FIOGC*(FJOTB*FNOTC-FJOTC*FNOTB) 

7  IF  (J-1)  A, 70, 69 

69  IF  (FOFI)  80,10,82 

80  SIGN2=-1,0 

GO  TO  83 

82  SIGN2=+1.0 

83  IF  (SIGN+SIGN2)  20,8,20 

70  IF  (FOFI)  6A, 10,66 
6A  SIGN=-1.0 

GO  TO  20 
66  SIGN=+1.0 

GO  TO  20 

8  DGG=1.0E-6 

IF  (DG-DGG)  10,10,9 

9  G=G+DG 
DG=DG/10.0 
GO  TO  20 

10  T=TT 
GG  =  G 

BB  =  SQRTF  (GG»»2  +  FK:(  I  )**2  ) 

RETURN 

END 
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C  FINDING  F»F.C.M»H 

SUBROUTINE  EFCMH ( L .A »B .C tNCVS* V 1 tX . Y ♦ ZZ .EE  » FF »CCK » DUMl . DUM2 »DUM3 ) 
dimension  X( 10) »Y( 10) »ZZ ( 10) ,EE( 10) .FF( 10) .CCK( lOtlO) »VI ( 10) 
DUM1»DUM2+DUM3 
DO  10  I=1»NCVS 
IF  ( I-l )  15.11»12 
15  PAUSE 
12  111*1-1 

DO  13  I  1=1.1  II 

IF  (VI ( I )-VI ( I  I  )  )  13.14.13 

13  CONTINUE 
II  T=X(I) 

G*Y( I ) 

BT=ZZ( I ) 

IF  (C)  1.2.3 

1  PAUSE  11111 

3  FN0TB=BESSELF(T*B.0.4) 

FN1TB*BESSELF( T*B.l .4) 

FJ0TC=BESSELF(T*C.0.1  ) 

FJ1TC=BESSELF(T*C.1.1 ) 

FN0TC»BESSELF( T*C.0.4 ) 

FN1TC*BESSELF(T*C.1.4) 

FIOGCsBESSELF(G*C.0.2) 

FI1GC*«ESSELF(G*C.1.2) 

2  FJ0TB=BESSELF( T*e.0.1 ) 

FJ1TB=BESSELF(T*B.1 .1 ) 

F10GA=BESSFLF(G*A.0.?) 

FI0GB»BESSELF(G*B.0.2) 

FK0GA=BESSELF(G*A.0.3 ) 

FK0G8=BESSELF(G*B.0.3) 

FI 1GB=BESSELF(G*B.1 .2 ) 

IF  (C)  4.5.6 

4  PAUSE  44444 

5  E=FJ0TB**2+FJ1TB**2 

FrFJ0TB**2«(  1.0+G**2/T**2+1.0/  (  (  T  ♦B*  (  F I  OGB«FKOGA-F  I0GA*FK.0G3  )  )**2  ) 
1) 

CK=2.0*( T*B*FI0GB*EJ1TB+G*B*FI 1GB»FJ0TB)/ ( (G**2*B«*2+T*#2»B*«2 ) *FI 
10GA*(Fj0Te**2+FJlTB**2) ) 

GO  TO  7 

6  BB=-( T*EJ1TC+G*EI 1GC«FJ0TC/F10GC) /(T«FN1TC+G»FI1GC#FN0TC/FI0GC ) 

E= ( B**2* ( ( E J0TB**2+FJ1TB**2 ) +2 • 0»BB» ( F J0Tb*FN0T3+F J1 TB*FNl TB ) +BB** 
12* (FN0TB**2+FN1 TB**2  » )-C**2*( ( F JOTC **2  +  F J 1 TC**2 ) +2 .0*BB* ( F JOTC*FNO 
2TC+FJ1TC*FN1TC ) +BB**2« ( FN0TC*»2+FN1 TC**2 ) ) ) / ( B**2-C**2 ) 

F=(  1 .0  +  G**2/T**2  )*  (B**2*(FJ0TB+BB*FiN0TB)**2-C**2*(  FJOTC+BB*FNOTC)* 
1*2)/(B**2-C**2)+(FJ0TB+BB*FN0TB)**2/(T**2*(B**2-C**2)*( FI0GB*FK0GA 
1-FI0GA*FK0GB)**2 ) 

CK=2.0*( (T*B*FI0GB*(EJ1TB+BB*FN1TB)+G*B*EI 1GB*(FJ0TB+BB*FN0TB) )-(T 
1*C*F IOGC*( FJl TC  +  BB*FN1TC )+G*C*F 1 1 GC* ( F JOTC+BB*FN OTC ) )  )  /  (E*FIOGA*(B 
2**2-C**2)*(G**2+T**2) ) 

7  EE(I)=E 
FF( I )=E 
CCK(L.I  )=CK 
GO  TO  10 

14  EE(I)=EE(in 
FF( I ) =FF  (  I  I) 

CCKiL.I  )=CCK(L.I  I  ) 

10  CONTINUE 
RETURN 
END 
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C  FREQUENCY  FUNCTION  SUBROUTINE 

SUBROUTINE  FREQ (Q»TU»DELO.NCVS  tFREAL, F IMAG ) 

DIMENSION  0(10) , TUI  10 ) ♦ FREAL ( 1 0 ) tFIMAGI 10) 

DO  10  I=1*NCVS 

CFCTR  =  0(1 )*(2.+DEL0+TU( I ) )*(DELO-TU( I ) )/ ( ( 1 .+DELO ) * ( 1 . +TU(  I  ) ) ) 
FREAL(I)  =  l./( l.+CFCTR**2) 

10  FIMAG(I)  =  -CFCTR»FREAL( I) 

RETURN 

END 


C  INTERPOLATION 

SUBROUTINE  TERPOL ( N, L ♦ Y .X . ANSR ) 
dimension  Y(3*10tlO) .ANSR(IO) 

DO  10  I=1*N 
A=Y(2.L, I ) 

B=(Y(3.L.I)-Y(l.L»n)/0,16 
C= ( Y( 1 ♦L.  I ) +Y( 3  tL.I )-2.0*Y(2*L ♦ 1) ) /0.0128 
ANSR( I ) =A+B»X+C*X**2 
10  CONTINUE 
RETURN 
END 


C  AGAIN  INTERPOLATION 

SUBROUTINE  A  INTER ( N, L * Y  ♦  X  ♦YY ) 

DIMENSION  Y(3*10.10) .YY(IO.IO) 

DO  10  1=1. N 
A  =  Y(2.L.  I ) 

B=(Y(3.L.I 1-Y( l.L.n )/0,16 
C= ( Y( l.L. I )+Y( 3.L.I )-2.0*Y(2.L. I ) )/0.0128 
YY(L.I)=A+B*X  +  C»X»«-2 
10  CONTINUE 
RETURN 
END 
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C  FORMATION  OF  MATRICES 

SUBROUTINE  FOMATRIN,  I  ,NCVS  ♦CICP  .CKM ,  FMP ,  FMM ,  FHP .  FHM  .  ZETAP  ,ZET  AM ,  YOG 
lOP.YOGOM.PPLUSR.PPLUSI .COLMR»COLMI .FLMPR.FLMPI ) 

DIMENSION  CKP( 10. 101 .CKM( 10 . 1 0 ) . FMP ( 10 , 10 ) . FMM ( 1 0 . 10 ) . FHP ( 10 . 1 0 ) . F 
1HM( 10.10) .ZETAPI 10.10) .ZETAMI 10.10) . YOqOP ( 10 . 10 ) .Y0G0M( 10.10) .PPLU 
2SR(20.20 ) .PPLUSI ( 20 .20 ) .COLMR I  20 . 20 ) .COLMI (20.20 ) . FLMPR ( 20 . 20 ) .FLM 
3PI (20.20) 

L  =  0 

NN=2»N 

DO  1  K=1  .NN.2 
J=1 
L  =  L  +  1 

PPLUSR( J.K) =CKP(L. I )*FMP(L.I )«FHP(L . I )*Y0G0P(L. I ) / ( 1 .0+ZET AP ( L . I ) ) 
PPLUSR( J.K+1 )=-CKM(L. I )*FMM(L. I )*FHM( L. I )*YOGOM(L. I ) / ( 1 . 0-ZETAM ( L . 
ID) 

PPLUSI (J.K)=0.0 
PPLUSI (J.K+1 )=0.0 

1  CONTINUE 
L*0 

DO  2  J=1.NN.2 
K  =  1 
L  =  L  +  1 

C0LMR(J,K)=CKP(L. I )*FMP(L.I )/(2.0*( 1 . O+ZE TAP ( L . I ) ) ) 

COLMR( J+1.K)=CKM( L.I )*FMM(L. I ) / ( 2 .0* ( 1 . O-ZETAM ( L . I ) ) ) 

COLMI (J,K)=0.0 
COLMI (J+1,K)=0.0 

2  CONTINUE 

DO  3  J=1 .NN 
DO  3  K=1.NN 
L  =  1 

FLMPR  (J.K)  =COLMR  (  J  .L  )  *PPLUSR  (  L  .  K  ) 

FLMPI  (  J.K)=COLMI  (  J.L)*PPLUSI  (L.O 

3  CONTINUE 
RETURN 
END 
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C  FORMATION  OF  D  MATR  IX 

SUBROUTINE  DMATRX (NCVSPl »N, 1 .BQPOS * BQMI N , BE ♦S^DR *0  I ♦ F IR ♦FI  I ) 
DIMENSION  BQPOSI 1 0 ♦  10 ) ♦BQMI N ( 10 ♦ 10 ) ♦BE ( 1 0 ) ♦ S ( 10 ) tDR ( 20 ^20 ) ♦ D I ( 20 ♦ 2 
10) ♦FIR(20^20) ♦FI  I (20^20) 

DO  1  J=1^N^2 
M  =  0 

DO  2  K=1 ♦N^2 
M  =  M+1 

IF  (J-K)  3^4^3 

3  DR(J^K)=0.0 
DI ( J^K ) =0.0 
DR( J+1 ♦< )=0.0 
DI ( J  +  1 ♦<  )=0.0 
DR( J^K+1 )=0.0 
DI ( J^K  +  1  )=0.0 
DR( J+1 ♦<+! ) =0.0 
DI  ( J  +  1 ♦<  +  ! ) =0.0 
FIR( J^K ) =0.0 

FI  I ( J.K ) =0.0 
FIR( J+1 ♦K)=0.0 
FI  I  (J+1 ♦K)=0.0 
FIR( J^K+1 ) =0.0 
FI  I  (  J.K  +  1 ) =0.0 
FIR( J  +  1 .K+l  )=0.0 
FI  I (J+1 .K+l  )=0.0 
GO  TO  2 

4  SXP=SINF (BE (NCVSPl )*S( I ) ) 

CXP=COSF (BE (NCVSPl )*S( I ) ) 

SYP=S I NF (BQPOS(M ♦NCVSPl )*S{ I ) ) 

CYP=COSF (BOPOS (M^NCVSPl )*S( I ) ) 

SYM=SINF (BQM IN (M ♦NCVSPl )*S( I ) ) 

CYM=COSF (BQM IN (M ♦NCVSPl )*S( I ) ) 

DR  ( J^K ) =SXP#SYP+CXP#CYP 

DI  ( J^K)=SYP*CXP-SXP#CYP 
DR ( J+1 ♦<+! ) =CXP*CYM-SXP*SYM 
13  DI ( J+1 ♦<+! )=-(SXP»CYM+SYM*CXP) 

FIR(J^K)=1.0 
FI  I ( J^K ) =0.0 
FIR (J  +  1 ♦<  +  !  )  =  1 .0 
FI  I (J+1 .K  +  l  )=0.0 
2  CONTINUE 
1  CONTINUE 
RETURN 
END 
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COMPLEX  MATRIX  MULTIPLICATION  (UP  TO  10X10) 

AXB  =  C 

N  =  ROWS  OF  A  AND  C 
M  =  ROWS  OF  B  AND  COLUMNS  OF  A 
L  =  COLUMNS  OF  B  AND  C 
SUBROUT INE  CMMP( AR,AI ♦N»M»BR*BI fL.CRtCI ) 

DIMENSION  AR( 20,20  )  ♦AI (20*20) ♦BRI 20*20 ) ,BI (20*20) .CR (20 »20) »CI (20» 
120 ) *D(20,20 ) ,E (20  *20  ) 

CALL  RMMP(AR*N *M*BR*L*D) 

CALL  RMMP(AI *N*M*BI *L*E) 

DO  10  I = 1  *  N 
DO  10  J=1*L 

10  CR(  I  » J) =D( I *J)-E( I *J) 

CALL  RMMP(AI *N*M*BR,L*D) 

CALL  RMMP(AR*N*M*BI *L*E) 

DO  20  I=1*N 
DO  20  J=1*L 

20  Cl ( l» J)=D(I *J)+E( I *J) 
return 
END 


C  COMPLEX  MATRIX  SUBTRACTION 

subroutine  CMSU3T(2GGL*FI *NN*SUM) 
dimension  ZGGL (20*20) *FI ( 20*20) *SUM(20*20) 
DO  1  J=1*NN 
DO  1  K=1*NN 

SUM(J»K)=FI ( J*K)-ZGGL(J*K) 

1  CONTINUE 
RETURN 
END 


C  REAL  MATRIX  MULTIPLICATION 

subroutine  RMMP  ( A*N*M*B*L*C) 
dimension  A(20  *20  )  *B( 20*20) *C(20*20 ) 

DO  5  1  =  1  *N 
DO  5  J=1,L 
(.(  I  *J  )  =  o.o 
DO  5  K=1*M 

5  C(  I*>^)  =  C(I*J)+A(  I  *K)*B(K,J) 

RETURN 

END 
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